Explorar separabilidad
if (!require("sf")) install.packages("sf")
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
if (!require("tidyverse")) install.packages("tidyverse")
## Loading required package: tidyverse
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.2 ✓ dplyr 1.0.6
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
if (!require("GGally")) install.packages("GGally")
## Loading required package: GGally
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
if (!require("np")) install.packages("np")
## Loading required package: np
## Nonparametric Kernel Methods for Mixed Datatypes (version 0.60-10)
## [vignette("np_faq",package="np") provides answers to frequently asked questions]
## [vignette("np",package="np") an overview]
## [vignette("entropy_np",package="np") an overview of entropy-based methods]
source("funiciones_sospechos.R")
e_zona <- st_read("Datos/poligonos_14_16.shp") %>% st_drop_geometry()
## Reading layer `poligonos_14_16' from data source `/home/alequech/Documents/usfsmex/capacitaciones/funcion_densidad_samof/Datos/poligonos_14_16.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1416 features and 59 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 2296366 ymin: 841959 xmax: 2340058 ymax: 872058.3
## Projected CRS: Lambert_Conformal_Conic
names(e_zona)
## [1] "fid" "unico" "OBJECTID" "OBJECTID_1" "mdmx_t1"
## [6] "mdmx_t2" "ipcc_dyna" "prob_cmb" "DNDVI" "MODO"
## [11] "SAMOF16" "Region" "Reg_Cuadr" "ID_UTE" "Val_T1"
## [16] "Val_T2" "Metodo" "Interprete" "OBS" "Shape_Leng"
## [21] "Shape_Area" "NDVI14" "NDVI16" "SUP_MNA" "Val_Cambio"
## [26] "oid_1" "change" "mdmx_t1_2" "mdmx_t2_2" "mdmx_t1t2"
## [31] "mdmx_cls" "ipcc_dyna_" "mdmx_t1_l" "prob_cmb_2" "maf_4"
## [36] "maf_5" "maf_6" "maf_1" "maf_2" "maf_3"
## [41] "tcd_t2" "tcd_t1" "met_t1_1" "met_t1_3" "met_t1_2"
## [46] "met_t1_5" "met_t1_4" "met_t1_7" "met_t1_6" "met_t2_6"
## [51] "met_t2_7" "met_t2_4" "met_t2_5" "met_t2_2" "met_t2_3"
## [56] "met_t2_1" "area" "id" "unico_raw"
# tcd tree cover density * 2015 en los casos posteriores
# mdmx_t1 etiqueta de cobertura t1
# mdmx_t2 t2
# ipcc_dyna_ etiqueta 1 a etiqueta 2
# prob_cmb_c prob. de cambio de t1 a t2 basado en la vecindad
# maf_5 coef. PCA
# met_t estadicas zonales de bandas
#t1 = 2014
#t2 = 2016
# Val_Cambio = C CAMBIO, p
my_colors <- c("P" = "#00b159","C" = "#d11141")
e_zona %>%
select(starts_with("met_") ,Val_Cambio) %>%
ggpairs(., aes(colour = Val_Cambio),
upper = list(continuous = wrap("points", alpha = 0.3)),
diag = list(discrete="barDiag", continuous = wrap("densityDiag", alpha=0.5 )),
lower = list(continuous = wrap("smooth", alpha = 0.3, size=0.1)),
legend = c(1,1))+
scale_color_manual(values = my_colors)+
scale_fill_manual(values = my_colors)
e_zona %>%
filter(SAMOF16 == 12)%>%
select(starts_with("met_") ,Val_Cambio) %>%
ggpairs(., aes(colour = Val_Cambio),
upper = list(continuous = wrap("points", alpha = 0.3)),
diag = list(discrete="barDiag", continuous = wrap("densityDiag", alpha=0.5 )),
lower = list(continuous = wrap("smooth", alpha = 0.3, size=0.1)),
legend = c(1,1))+
scale_color_manual(values = my_colors)+
scale_fill_manual(values = my_colors)
Función de densidad - 1 paso
Todas las categorias
variables <- names(e_zona)[startsWith(names(e_zona),"met_")]
formula <- formula(paste("~", paste(variables, collapse = "+")))
bw <- np::npudensbw(formula, data = e_zona, bwmethod="normal-reference")
f <- npudens(bw)
df_sospechosos <- cbind(e_zona, f$dens)
df_sospechosos$sospechoso <- with(e_zona, ifelse(f$dens <= quantile(f$dens, 0.10), 1, 0))
ct_t <- table(df_sospechosos$Val_Cambio, df_sospechosos$sospechoso)
ct_t
##
## 0 1
## C 108 6
## P 1165 137
prop.table(ct_t) # cell percentages
##
## 0 1
## C 0.076271186 0.004237288
## P 0.822740113 0.096751412
prop.table(ct_t, 1) # row percentages
##
## 0 1
## C 0.94736842 0.05263158
## P 0.89477727 0.10522273
prop.table(ct_t, 2) # column percentages
##
## 0 1
## C 0.08483896 0.04195804
## P 0.91516104 0.95804196
ct <- table(df_sospechosos$Val_Cambio,df_sospechosos$sospechoso)
ct
##
## 0 1
## C 108 6
## P 1165 137
prop.table(ct) # cell percentages
##
## 0 1
## C 0.076271186 0.004237288
## P 0.822740113 0.096751412
prop.table(ct, 1) # row percentages
##
## 0 1
## C 0.94736842 0.05263158
## P 0.89477727 0.10522273
prop.table(ct, 2) # column percentages
##
## 0 1
## C 0.08483896 0.04195804
## P 0.91516104 0.95804196
para una categoria
sb <- e_zona %>%
filter(SAMOF16 == 12)
bw_sb <- np::npudensbw(formula, data = sb, bwmethod="normal-reference")
f_sb <- npudens(bw_sb)
df_sospechosos_sb <- cbind(sb, f_sb$dens)
df_sospechosos_sb$sospechoso <- with(sb, ifelse(f_sb$dens <= quantile(f_sb$dens, 0.10), 1, 0))
ct_sd <- table(df_sospechosos_sb$Val_Cambio,df_sospechosos_sb$sospechoso)
ct_sd
##
## 0 1
## C 18 7
## P 315 30
prop.table(ct_sd, 1)
##
## 0 1
## C 0.72000000 0.28000000
## P 0.91304348 0.08695652
prop.table(ct_sd, 2)
##
## 0 1
## C 0.05405405 0.18918919
## P 0.94594595 0.81081081
# densidad_1paso <-
# function(df_cat,
# poutlayer = 0.10,
# var = c("navevar1", "navevar2")) {
# formula <- formula(paste("~", paste(var, collapse = "+")))
# bw <- np::npudensbw(formula, data = df_cat[,var], bwmethod="normal-reference")
# #f <- npudens(tdat = df_cat[, var])
# f <- np::npudens(bw)
# df_cat$fdens <- f$dens
# df_cat$sospechoso <- with(df_cat, ifelse(df_cat$fdens <= quantile(df_cat$fdens, poutlayer), 1, 0))
# return(df_cat)
# }
var_t1 <- names(sb)[startsWith(names(sb),"met_t1")]
sb_t1 <- densidad_1paso(sb ,poutlayer = 0.10, var = var_t1)
ct_sd_1 <- table(sb_t1$Val_Cambio,sb_t1$sospechoso)
ct_sd_1
##
## 0 1
## C 20 5
## P 313 32
prop.table(ct_sd_1, 1)
##
## 0 1
## C 0.80000000 0.20000000
## P 0.90724638 0.09275362
prop.table(ct_sd_1, 2)
##
## 0 1
## C 0.06006006 0.13513514
## P 0.93993994 0.86486486